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<N ■ Abstract 

> . 

We present numerical results for the phase diagram of lattice QCD at finite tem- 
perature in the formulation with twisted mass Wilson fermions and a tree-level 
CO ' Symanzik-improved gauge action. Our simulations are performed on lattices with 

I/-) ■ temporal extent N T = 8, and lattice coupling f3 ranging from strong coupling to 

the scaling domain. Covering a wide range in the space spanned by the lattice 
coupling (3 and the hopping and twisted mass parameters k and fi, respectively, we 
obtain a comprehensive picture of the rich phase structure of the lattice theory. In 
particular, we verify the existence of an Aoki phase in the strong coupling region 
and the realisation of the Sharpe-Singleton scenario at intermediate couplings. 
■ In the weak coupling region we identify the phase boundary for the physical fi- 

nite temperature phase transition/crossover. Its shape in the three-dimensional 
parameter space is consistent with Creutz's conjecture of a cone-shaped thermal 
transition surface. 



1 Introduction 



The quark-hadron transition predicted by Quantum Chromodynamics (QCD) is a central 
aspect in describing strongly interacting matter at finite temperature T. The early 
universe has passed through this transition, which is now hoped to be reproduced in 
ultrarelativistic heavy ion collision experiments. With a transition temperature T c ~ 
150 — 200 MeV, the relevant energy scale is in the inherently non-perturbative regime of 
QCD, calling for numerical simulations of lattice QCD. We consider net baryon density 
zero, for which Monte Carlo methods are directly applicable (cf. [H, 0] for a review). 
Most previous investigations of lattice thermodynamics have used the Kogut-Susskind 
(KS) staggered fermion formulation or the Wilson fermion formulation, amended by the 
improvement programme originally laid out by Symanzik 0, 0, [i], @| . Both approaches 
have led to predictions for the nature of the transition, the transition temperature and 
the equation of state 0, S, 0, 0, 11, 12|. In particular, based on simulations of a 
staggered action ref. [lOj states that the finite temperature transition in continuum 
QCD at the physical point is analytic with a transition temperature depending on the 
observable it is extracted from. Given these results on the one hand and the continuing 
controversy on the universality of the staggered formulation 13|, ll4| on the other, a check 
of these findings with fermions of the Wilson type appears desirable - even though the 
phase structure of the latter is more complex compared to KS fermions due to explicit 
chiral symmetry breaking at 0(a). 

An economic alternative to the computationally involved Symanzik improvement 
and/or smearing techniques for the Wilson fermion formulation might be Wilson twisted 
mass (Wtm) fermions 15j,[l6|. This formulation allows for automatic (9(a)— improvement 
of physically relevant operators by merely tuning Wilson's hopping parameter k to its 



critical value (so called maximal twist, cf. [17| for a detailed review), without computa- 
tional overhead. This fermion discretisation has been studied extensively at zero tem- 
perature by the European Twisted Mass Collaboration (ETMC) \l~ 



IS, 20, 21, 22, 23 



providing a growing amount of data relevant for tuning to maximal twist, renormalisation 
factors and setting the scale 2J, |25j, |26|, |27| . In particular, simulations at pseudoscalar 



masses mps ~ 300 MeV are now feasible, and considerable progress in controlling cut-off 
and finite size effects has been achieved by using lattice chiral perturbation theory [28 . 
The price to pay for automatic improvement is a more complex phase structure of the 
lattice model due to the introduction of an additional parameter, the twisted mass /io, 
as well as explicit flavour symmetry breaking by the twisted mass term for /io ^ 0. 

The aim of our present study is to investigate the applicability of Wtm fermions 
for lattice thermodynamics with two degenerate quark flavours. For intermediate quark 
masses only a smooth crossover is expected in this case, which should turn into either 
a second or first order phase transition as the chiral limit is approached. The latter 
situation is, of course, very difficult to establish and currently still under debate, for 
some of the most recent attempts with opposing conclusions, see 29|, |30j . In this work, 



we restrict ourselves to a study of the phase structure of the Wtm fermion model with 
a temporal lattice size N T = 8. This is a prerequisite to establish whether simulations 
in the physically interesting regime and using automatic 0(a) improvement are feasible 
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without running into unphysical phases. The results we present here answer this question 
in the affirmative. 

In Sec. [2] we specify the action and observables used in our simulations. The current 
knowledge about the phase structure for Wilson type fermions at finite temperature and 
that of Wtm fermions at zero temperature is summarised in Sec. [3j We present our 
simulation results in Sec. H] and discuss them in the light of lattice chiral perturbation 
theory. Finally, Sec. [5] contains a discussion and our conclusions. Preliminary results of 



this work have already been reported in [3l|, |32j, [33 



2 Action and observables 

Our setup consists of a hypercubic lattice N% x N T with periodic boundary conditions 
for the gauge fields, whereas the fermionic fields are subject to (anti-) periodic boundary 
conditions in the (time) space direction. We use the tree-level Symanzik improved gauge 
action given by 



S g [U]=P 



c ° E i 1 - ^Re(Tr[tf(p)])) + a J2 (i - lMnu(R)])) 



(1) 



where (3 = with g the bare gauge coupling and U(P),U(R) denote the path- 

ordered products of link variables along closed loops of length 4a (plaquettes) and 6a 
(rectangles), respectively. The weight coefficients satisfy cq + 8c\ = 1 (normalisation) 
and Ci = —1/12 (tree- level improvement condition). This choice of gauge action is 
motivated by allowing us to exploit ETMC's data for T = 0. 

The Wilson twisted mass action for the fermion sector in the so-called twisted basis 
reads 

S F [t/j , $ , U] = a 4 (Dw + m + i/io7sr 3 ) ip(x)] , (2) 

with D w = 7^2^ (y v + ^*) ~~ t^m^m the Wilson-Dirac operator. This action can be 
rewritten as 

S F [V , * , U] = a 4 (7, Ya (V, + V^) - ^e"^- 3 V^V, + M )*(i)] (3) 

X 

in the so-called physical basis ^, which is related to the twisted basis z/>, ip by the 
non-anomalous chiral rotation 



\P = exp (^757" 3 ) i> 



^ = ^exp ^| 75 r 3 ) , (4) 

with 

M = \Jml + nl , tan(w) = fi /m (5) 
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denoting the bare polar quark mass and the bare twist angle, respectively. Finally, 
introducing Wilson's hopping parameter k = l/(2am + 8r) and rescaling the fermion 
fields according to ip — > a/ a 3 / (2K)ip leads to the familiar form 

S f [*Ij ,*fj,U] = J2 [^0) (! + i^a^r 3 ) ip(x) 

(6) 

- « ^ if)(x) ((r - 7 M )f/ M (x)^(x + p) + (r + l^)U^(x - afi) ] i)(x - ft)) 



Because of the spin-flavour structure of the twisted mass term, ia / u 7 5 r 3 , parity is only 
a symmetry when combined with a discrete flavour rotation or a sign change /xq — * —^o, 
while flavour symmetry is broken explicitly according to the pattern S77y(2) — > C/(3)(l), 
where C/(3)(l) is the subgroup of flavour rotations generated by r 3 . 

For our numerical simulations we use a generalised hybrid Monte-Carlo algorithm 
with even-odd preconditioning 34 




and multiple time scale 
(cf. [3^] for algorithmic 



Hasenbusch trick 
integration according to the Sexton- Weingarten scheme |3j 
benchmarks and 39j for a more recent collection of algorithmic improvements). 

In our investigation of the phase structure we measure and analyse the following 
observables: 



1. (averaged) plaquette P 

1 



P 



6N*N T 



(7) 



X /!<!/ 



2. (real and imaginary part of the) Polyakov loop L 

'N T -1 



— VTr 



3iV 3 



n 

i=0 



3. their susceptibilities according to 

X (0) = iV 3 ((O 2 ) - (O) 2 ) , (9) 

4. fermionic observables: scalar condensate (ipip), the order parameter for parity- 
flavour symmetry breaking (ipi'-f 5 T 3 tfj) and the pion norm 

|tt| 2 = ^(^( x ) 75 ^^(a;)^(0)75 > ( 10 ) 

X 

where = r 1 ± it 2 . 

Discontinuous phase transitions are in principle indicated by any observable, but depend- 
ing on the particular dynamics one may be more sensitive than the other. In particular, 
the Polyakov loop is an order parameter for Z(N) -breaking and hence susceptible to 
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Figure 1: Left: Aoki phase and thermal transition K t (f3) for untwisted Wilson fermions. 
/Scusp denotes the tip of the Aoki phase and /3 qu the critical coupling for the deconfinement 
transition in the quenched limit. Right: Aoki phase and first order bulk transition for 
twisted mass Wilson fermions at zero temperature [17( . 



signals of deconfinement, which get more pronounced the heavier the quark mass. On 
the other hand, the pion norm and scalar condensate are sensitive to small eigenval- 
ues of the Wtm operator and hence detect signals of chiral symmetry breaking in the 
light mass regime. The pseudo-scalar condensate is the order parameter for spontaneous 
parity-flavour symmetry breaking in the plane {/i = 0} of the (/?, k, /x )— parameter 



space 



40, 41, 42 



3 Current knowledge of the phase structure 



3.1 Wilson fermions at T ^ 

The phase structure for standard and improved Wilson fermions is determined by two 
phenomena: spontaneous parity-flavour symmetry breaking and the finite temperature 
transition. The existence of a phase of spontaneously broken parity-flavour symmetry, 
commonly called Aoki phase, was introduced theoretically and verified numerically in 
(iol lill HI, 0, Hl| . Refs. [H, 0, 0, Hi] demonstrated that the Aoki phase is located 
in the strong coupling region and ends in a cusp at a finite ft cusp , from which the chiral 
critical line k c (/3) emanates and continues towards the weak coupling limit, cf. Fig. [T] 
(left). The Aoki phase is repeated for larger values of k, corresponding to the fermion 
doubling phenomenon. The value of /3 cusp depends on the choice of action and the 
temporal lattice extent. Being a lattice artefact, the Aoki phase does not affect physics 
in the weak coupling region. 

The finite temperature transition or crossover is characterised by a physical transition 
temperature which depends on the renormalised quark mass, T c (m#), and therefore, 
through rriR, on the bare quark mass. Since T = l/(a(/3)N T ) on the lattice and the bare 
quark mass is parametrised in terms of k, pairs of (rriR, T c ) translate to pairs (n t (Pc), Pc) 
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on the lattice describing the phase boundary for finite temperature transitions at fixed 
N T . This line has been mapped out in detail by the CP-PACS 0, H Q and DIK 
@, IHi 52] collaborations for the case of an (9(a)-Symanzik-improved fermion action 
on one hand and the standard plaquette or the renormalisation group improved gauge 
action on the other with the following results: 

1. The thermal line for fixed N T runs from a finite /3 qu in the quenched limit (n t = 0), 
where the transition is of first order, towards the zero temperature chiral critical 
line K c (f3,T = 0) and lower (3. This is consistent with the physical expectation 
that the transition temperature increases with the quark mass. 

2. In the vicinity of k c ({3,T = 0) and the cusp of the Aoki phase, the thermal line 
appears to run close and parallel to the critical line; despite large computational 
efforts a clear separation of the thermal line from the border line of the Aoki phase 
has not yet been achieved on N T = 4 lattices 48| in case of the pure Wilson gauge 
action. 



3.2 Twisted mass fermions at T = 

The zero temperature phase structure of Wilson twisted mass fermions has been inves 



19, Ml 21 



23 



tigated in a series of publications by the ETM Collaboration 
The parameter space is extended from the {/3,k}— plane to the 3d {/3, k, fio}— space. 
Hence, in the strong coupling region an Aoki phase should exist in the (3 — k— plane for 
< (3 < /5 C usp, cf. Fig. [1] (right). Adjacent to it a surface of first order phase transi- 
tions extends into the weak coupling region. It is perpendicular to the Aoki phase and 
includes the chiral transition line at zero twist. This surface is predicted from chiral 
perturbation theory 53|, |54| and has been established in zero temperature simulations 



[18( | . This transition for ^ > is also a lattice artefact, correspondingly its width in 
the fiQ— direction decreases with increasing f3 and can be made smaller at a fixed value 



of (3 by using improved gauge actions 20]. Finite temperatures correspond to finite N T 
on the lattice, which will amount to small shifts of the transition surface compared to 
T = 0. 

These unphysical first order transitions and their suppression are a crucial issue for 
simulations at small quark and corresponding pseudoscalar masses. At maximal twist 



25j, |26j, realised by setting k = K c (f3), the width of the surface in /io-direction implies a 
lower bound /xo c (/5) for the twisted mass and hence for the quark mass (cf. Eq. ([5])) at 
fixed j3, in order to avoid unphysical phase transitions. 



3.3 Twisted mass fermions at T ^ 0: Creutz's cone conjecture 

In F 



55, 56 



Creutz conjectured a finite temperature phase structure for untwisted and 
twisted Wilson fermions based on continuum symmetry arguments. 

In the continuum theory the twist angle tan(u;) = fto/mo is an irrelevant parameter 
and the rotation exp(zu;/275T 3 ) a mere change of integration variables in the path inte- 
gral, leaving the path integral measure invariant. Hence, T c can only be a function of 
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Figure 2: Schematic view of Creutz's cone conjecture, derived from [56 . 



the renormalised quark mass Mr = ym R + fi R , i. e. T c is constant on lines of constant 
Mr, which correspond to circles. In bare parameter space, because of the renormali- 
sation factors for the masses these curves are ellipses. This picture should translate to 
the lattice, modified by cut-off effects. For fixed (3, the lines of constant Mr are closed 
curves that become ellipses in the continuum limit. Hence, for fixed N T there should be 
an interval /? CUS p < (3 < (3qu where the surface of finite T transitions wraps around the 
chiral transition line. Since T c grows with Mr, this surface should have a conical shape 
opening in the direction of growing (3. Fig. [2] gives a schematic summary of the expected 
phase structure. 

While the existence of a finite temperature transition for every value of the twist 
angle is dictated by the continuum theory, the strength and nature of the transition 
are subject to discretisation effects and thus generically dependent on the twist angle 
uj. This will lead to a distortion of the ellipse and conical shape. For sufficiently light 
quarks, this expectation can be made more quantitative in the framework of lattice chiral 
perturbation theory. At next-to-leading order (NLO), the expression for the pion mass 



is 



57, 58 



7T- 1 - 



M> + j 2 ((2L e 



- L 45 ) {M'y + M'a cos{u) (2W - W) + 2a 2 cos 2 {u)W) (11 



(M') 2 ( M' 



where the quark masses are now renormalised ones, 



fi = Z^[i , m = Z, 
2B Z^ , rh' = 2BqZ, 



m - 

m (m 



Furthermore, M' = \J (i 2 + (m 1 ) 2 , a = 2W a, and \x 
The 5's, Ws, L's can in principle all be fixed by fitting lattice chiral perturbation theory 
to zero temperature simulation results. Currently, this is not completed yet. Once the 
NLO formula describes the data sufficiently well, knowledge of one thermal transition 
point on the cone will suffice to predict the entire distorted ellipse for a given (3 and pion 
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Figure 3: Fisher plots for /? = 1.8 (upper left panel), (3 = 3.0 (upper right panel) and 
(3 = 3.4 (lower panels) for k < k c ps 0.1825 (left) and k > n c (right). 

mass. 

4 Simulation results 

The goal of our numerical simulations is to establish and locate the different phases 
discussed in the previous section. The following analysis of our numerical data thus 
proceeds along the (3— axis from strong coupling (/? < 1.8) to the scaling region {(3 > 3.8). 
The parameter sets for which simulations have been carried out are quoted in Tab. [T]in 
the Appendix. Throughout this paper we have used lattices with temporal size N T = 8 
and spatial size N c = 16, except in a few cases where it will be indicated explicitly. 

4.1 Strong coupling: the Aoki phase 

We checked for the existence of an Aoki phase at (3 — 1.8, 3.0, 3.4. For the analysis 
we consider so-called Fisher plots showing (ipi^T^ip) 2 as a function of 2k/i / ' (xjjijsT 3 ^) 
(cf. (42! and references therein for a motivation of this approach). One then inspects the 
behaviour in the limit 2/t/i / (ipi'-f^T 3 '^) — > 0: for (j3, k) points outside the Aoki phase 
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K K 

Figure 4: Metastability signals in the order parameter (ipi^sT 3 ^) (left) and the Polyakov 
loop (right) in the neighbourhood of the critical line, /x ~ 0.0068. The lines were added 
for visual guidance. 



the curves will hit the abscissae at a non- vanishing value of 2k^o/ (^75 r3 ^); a point on 
the phase boundary would lead to a straight line through the origin and a point within 
the Aoki phase would show a positive intercept with the ordinate axis. Our results are 
shown in Fig. [3J 

The upper left panel shows the situation for the strongest coupling we considered, 
(3 = 1.8. We checked for finite-volume effects on lattices with N a = 16, 24, 32 and found 
them to be negligible within the statistical uncertainty. From the curves we conclude 
that the point (J3, k) = (1.8, 0.24) lies inside the Aoki phase, whereas ((3, k) = (1.8, 0.23) 
is outside. At (3 — 3.0 the upper right panel suggests that the Aoki phase has narrowed 
to the interval 0.203 < n < 0.204. The curve for k = 0.2035 was measured for N a = 16 
only. The point (j3, hi) = (3.0, 0.2035) appears to be a candidate for spontaneous parity- 
flavour symmetry breaking. 



4.2 Intermediate coupling: the Sharpe-Singleton scenario 

The situation changes notably when proceeding to (3 = 3.4, as shown in the lower panels 
of Fig. [3j Neither below nor above n c ~ 0.1825 do we see any indication of spontaneous 
symmetry breaking. Hence, within our resolution the existence of an Aoki phase is ruled 
out for (3 = 3.4. This assertion is supported by the emergence of another phenomenon. 
We find clear signals of metastability in all observables when scanning in k across the 
region of the chiral critical line, provided /io was chosen small enough. Fig. H] shows 
this behaviour for two different observables. We interpret these two-state signals as tied 



to the surface of first order phase transitions discussed in Sec. 13.21 (Note that in [59 
metastability signals were also reported in simulations with pure Wilson fermions at 
finite temperature and iV r = 6). 
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Figure 5: Real part of the Polaykov loop (left) and the pion norm (right) around k c at 
P = 3.65 and /x = 0.0068. 
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Figure 6: Integrated autocorrelation time for Re(L) as a function of k for j3 = 3.65, /io 
0.0068 (left) and (3 = 3.75, /x = 0.0070 (right). 



4.3 Weak coupling: finite temperature transition near k c 

In this section we discuss our observations of the finite temperature transition, specifi- 
cally its evolution along the j3— axis. Despite the reasonable statistics of typically C(15K) 
trajectories per data point, our results are still noisy and remain somewhat qualitative. 
This is also due to the fact that we are, most likely, dealing with a very soft crossover 
rather than a true phase transition. 

Let us start by reconsidering the right panel of Fig. HI The real part of the Polyakov 
loop rises significantly when k is increased, indicating the transition to a deconfined 
regime. Upon crossing k c , however, it enters the metastability region discussed in the 
previous section and then drops again. We interpret the rise and fall of the Polyakov loop 
along the qualitative expectations from Fig. [2] as an indication for two finite temperature 
transitions: first a confinement — > deconfinement transition when approaching n c from 
below, followed by a deconfinement — ► confinement transition at k > K c . This finite 
temperature transition is masked by the unphysical bulk transition around n c if /i is 
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Figure 7: The real part of the Polyakov loop (left) and its susceptibility (right) for 
(3 = 3.75, 3.775, 3.8, 3.9 and /i = 0.005. The vertical lines mark n c {T = 0,j3) for 
P = 3.9, 3.8, 3.75 from left to right. 



chosen too small. However, according to Sec. I3.2I this should be disentangled at larger 
values of j3. 

Indeed, choosing (3 = 3.65 while keeping /i constant, we observe the behaviour 
displayed in Fig. [51 The metastabilities are now suppressed and (Re(L)) peaks before 
gently declining again. Similar signals are obtained from other observables, with a 
sharper peak in the pion norm marking the chiral transition at k c , a remnant from the 
zero temperature bulk transition. 

An inspection of the data points and their statistical uncertainties reveals an asym- 
metry between k < n c and k > k c . This is because for k > k c the spectrum of the 
twisted mass Dirac operator is shifted in the negative real direction, which implies the 
occurence of very small eigenvalues and negative quark masses The results of an 



autocorrelation analysis are shown in Fig. The integrated autocorrelation time T mt 
of the Polyakov loop rises notably when passing k c , and it continues on a heightened 
plateau at negative quark masses. Hence, while we believe to have good evidence for its 
existence, the upper half cone of thermal transitions n t > n c is not easily accessible by 
simulations. 

By a further increase of (3 to (3 > 3.75 we approach the physically interesting weak 
bare coupling region. Looking at the left panel of Fig. [71 we once more find the real 
part of the Polyakov loop displaying a maximum around K c , which is consistent with 
the cone structure. Note also that, in agreement with qualitative expectations, the 
cone is widening at larger (3. It is customary to locate phase transitions more precisely 
by generalised susceptibilities, which mark the points of maximal fluctuations. These 
are shown on the right hand side of Fig. [7] for the Polyakov loop. The first transition 
appears to weaken with decreasing j3 and turns into a rather flat shoulder at (3 = 3.75. 
Decreasing (3 implies decreasing T c , and since T c is an increasing function of quark 
mass, the transitions for decreasing /3-values belong to a sequence of decreasing mass 
parameters. Our signals are thus consistent with a weakening of the deconfinement 
transition when moving away from the quenched regime. The latter was roughly located 
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Figure 8: Plaquette susceptibility (left) and pionnorm (right) for (3 = 3.75, 3.775, 3.8 
and /i = 0.005. 



for the tree-level Symanzik improved gauge action by p erforming a series of (3— scans at 
very large quark mass and /io = 0.005. As expected 45J the finite temperature transition 
line bends down from the chiral chritical line towards the j3— axis and hits the latter at 
a finite value (3 qu > 4.50(5) (cf. Table ED- 

Fig. [8] shows that the qualitative picture obtained from the Polyakov loop is repeated 
in the plaquette susceptibility (left) and the pion norm (right). Note that in all cases 
the location of the first thermal transition is at n t < n c , consistent with the qualitative 
expectations in Figs. [TJ [2J This becomes more obvious in Fig. O where we compare in 
detail the relative position of k c and n t for the confinement-deconfinement transition for 
(3 = 3.75, 3.8, as observed in all three observables we consider. Whenever a peak in a 
susceptibility can be isolated, the value for K t can be consistently determined from any 
observable by fitting Gaussian curves through the peaks. Our estimated values for the 
critical couplings are collected in Table 

An additional run on a larger lattice does not seem to result in a growth of the 
susceptibility peaks, as displayed in Fig. [101 a behaviour consistent with a crossover. 
Our more detailed study at (3 = 3.75, which will be reported in Sec. I4.51 will indeed 
indicate a pion mass over 400 MeV, for which only a smooth crossover behaviour is 
expected. 

The situation changes drastically for k > k c : the persistence of the Polyakov loop 
susceptibility at a high level — a phenomenon we again relate to the peculiar dynamics 
and its algorithmic problems in this region of quark masses — leaves the localisation 
of the second transition to be a delicate task. However, we consider the descent of the 
Polyakov loop as a clear signal for the deconfinement-confinement transition at k > k c . 
We further tested this interpretation by a (3— scan with k = 0.17 above k c . A signal for 
a transition was indeed detected and is also given in Table [2J 

In contrast to the situation at k < k c , these signals are accompanied by a signif- 
icantly weaker or altogether missing corresponding signal in other observables. This 
demonstrates the influence of cut-off effects in the simulation of the lattice theory: while 
the thermal transition in the continuum is manifestly invariant under changes of the 
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twist angle, we find this invariance broken by lattice artefacts such that positive and 
negative values for the twist angle can result in the two rather distinct signal scenarios 
we encountered. Let us recall, however, that the physically interesting region is the lower 
half of the cone with k < k c , where the signal for the transition can be located and is 
consistent among different observables. 



4.4 Second branch of finite temperature transitions at large k 

In order to gain a global perspective of the phase structure, we extended our investigation 
towards k ^> k c , even though this region is not relevant for continuum physics. The 
results for the Polyakov loop and its susceptibility of several k— scans are collected in 
Fig. [TTJ There appears to be another line of potential confinement — > deconfinement 
transitions moving towards the chiral critical line with increasing (3. The existence of 
such a branch of transitions is consistent with the conjecture in Fig. [2j it matches the 
line bending upwards from the upper end of the cone at large (3 towards the first doubler 
region. In order to give a rough idea of the distance between the two branches we 
present the data over a large K-range including also the data near k c for (3 = 3.75. The 
distributions of the Polyakov loop in Fig. [12] and especially the pronounced metastability 
found when simulating on a spatially enlarged lattice with hot and cold starts suggest a 
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Figure 11: Real part of the Polyakov loop (left) and its susceptibility (right) for (3 = 
3.4, 3.45, 3.65 with /i = 0.0068 and (3 = 3.75 with /x = 0.005. The vertical line marks 
K c (f3 = 3.75). Lines on the right are added for visual guidance. 
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Figure 12: Distributions of Re(L) for (3 = 3.6, /xo = 0, k — 0.22. 
possible first order nature of the transition in this parameter regime. 
4.5 Thermal transition line at fixed (3 

Having established the overall qualitative phase structure, let us now focus on the shape 
of the thermal transition surface in a plane of fixed f3. To this end we supplement 
our K-scan at (3 = 3.75, /x = 0.005 from Figs. CD, El by additional ones at larger /x 
values. The qualitative behaviour remains unchanged, but the two thermal transitions 
should approach each other with increasing /i , as expected for a slice of a conical 
shape. Eventually, for large values of /xo one is outside the cone and the transition signal 
should be lost. This is shown in Fig. [T3] (left), where the two transitions can be seen in 
the rise and fall of the real part of the Polyakov loop. For the largest twisted masses 
Ho = 0.025, 0.035 this behaviour is definitely lost. 

In order to locate the transition points more precisely, we have attempted to fit 



13 



0.013 
0.012 
0.011 
0.01 
0.009 
0.008 
0.007 
0.006 
0.005 
0.004 
0.003 



fi =0.065 
M-n=0.010 
M- =0.014 
fi =0.025 
fi =0.035 



ft 



1 



0.011 
0.01 
0.009 
0.008 
0.007 
0.006 
0.005 
0.004 



LO formula 
NLO formula? 
" %(Re(L)) 

Sinf 



0.165 



0.166 



0.167 



0.168 



0.165 0.1655 0.166 0.1665 0.167 0.1675 0.168 

K 



Figure 13: Left: Polyakov loop expectation value for different values of /j,q. For /x > 
0.025, the transitions are lost. Right: Data points represent the phase boundary at 
(3 = 3.75, estimated from the peaks of susceptibilities. The leading order part of Eq. (fTTjl 
does not fit the data, the NLO-formula contains undetermined constants that have been 
estimated by order of magnitude only. The vertical line marks k c (/3 = 3.75), i.e. maximal 
twist . 



Gaussian curves through the peaks in the susceptibilities, as shown in the examples in 
Fig. O In particular, we may use the clear first transition in Fig. Q3] (left) for a rough 
estimate of the corresponding pion mass. According to Table [2] the thermal transition 
point Kt is in this case very close to maximal twist K c (j3 = 3.75), for which llq = 0.005 



corresponds to a pion mass ~ 400 MeV [61J. Since smaller n t implies a slightly larger 



mass, the thermal cone corresponds to a pion > 400 MeV. 

As was discussed before, different observables display signals of different quality, and 
our data are not precise enough to do an unambiguous determination of both transitions 
in all observables. Moreover, the underlying physical transition is a broad crossover. 
For parameter values where entry and exit of the deconfined phase are close, the signals 
overlap and Gaussian fits do not work reliably, cf. Fig. [T41 (right). Nevertheless, in those 
cases that could be resolved, the peak positions are consistent between all observables, 
so we have used the observable with the strongest signal in each case. This results in a 
qualitative picture of the phase boundaries as shown in Fig. [13] (right). 

As discussed in Sec. I3.3[ in principle the cone of thermal transitions corresponding 
to a given pion mass at fixed j3 can be described by lattice chiral perturbation theory, 
Eq. (1111) . provided the mass and cut-off effects are sufficiently small. We observe that 
the leading order formula with only the first term in Eq. ([711) is unable to accomodate 
our data, even if an unknown constant is left open as a fitting parameter. This is not 
surprising given the large pion mass corresponding to the transition at (3 — 3.75. For the 
NLO curve in Fig. [13] (right) we employed some very rough estimates for the unknown 
parameters in Eq. (11 II) to check for consistency with that formula. It would now be 
interesting to really determine the unknown constants in the full NLO-formula Eq. (1111) . 
in order to see whether it correctly describes the data. 
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Figure 14: Pion norm (left) and Polyakov loop susceptibility (right) showing two sub- 
sequent transitions at j3 — 3.75, /io = 0.005. 

5 Discussion and conclusions 

We investigated the phase structure of two-flavour Wilson twisted mass fermions at 
non-zero temperature with special interest in the thermal transition close to the zero 
temperature chiral critical line. Our estimates for transition points from numerical 
simulations are summarised in Fig. [T3J, augmented by data for k c (/3 > 3.75) from the 



ETM collaboration [24|, |61| and the results of a rough tracking of the thermal line towards 
the quenched limit. 

In Fig. [16l we present an expanded view of the scaling region. The values of r^T 



are shown along the upper horizontal axis (cf. reference [2ljj for the value of r /a in the 
chiral limit). Note that for improved staggered and Wilson fermions, extrapolations of 
the transition temperature to the chiral limit give tqT c ~ 0.45 ji, 62], much smaller than 
our explored range. This is consistent with our large value of the pion mass and shows 
that larger N t > 10 are required to approach the chiral limit of two flavor QCD. 

Using the tree-level Symanzik improved Wilson gauge action, we produced substan- 
tial results in support of both the existence of an Aoki phase at strong coupling and the 
realisation of the Sharpe-Singleton scenario adjacent to the Aoki phase, with a bulk tran- 
sition surface narrowing towards the weak coupling region. Our results for the thermal 



transition are consistent with the conical structure proposed in [56( . 

From our extensive studies exploring a wide range of parameter space we may then 
conclude that the global phase structure of twisted mass lattice QCD is understood. 
Moreover, we have shown that for lattices with N T > 8 this phase structure is beginning 
to show the features expected in the continuum. In particular, it is then feasible to tune 
to maximal twist by setting k = k c {(3), thus ensuring automatic 0(a)-improvement, and 
investigate the thermal transition in the {/?, /io}-plane without running into unphysical 
phases. A recent perturbative study has shown that in this case maximally twisted 
Wilson fermions show comparable scaling behaviour in the high temperature region as 



standard staggered fermions |63j, while such simulations should not be any more compli- 



cated than those with untwisted Wilson fermions. It would thus be interesting to apply 
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Figure 15: Summary of numerically obtained transition points for < /j,q < 0.007. The 



data for k c (/3), (3 > 3.75 are taken from 24|, |61 
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Figure 16: The region (3 > 3.75. The lines are splines to guide the eye where the 
quenched and zero coupling limits have been taken into account. 
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the twisted mass formulation to systematic investigations of QCD thermodynamics. 

With results in this direction we hope to come up soon. Employing the twisted mass 
formulation we would like to add an independent universality check to the resolution 
of the ongoing controversary concerning the value of the critical temperature(s) and 
the issue of a possible separation between the deconfming transition and the chiral 
transition (for a recent paper on this subject see [64J and papers cited therein). We 
are confident that the automatic 0(a) improvement will be technically advantageous for 
this task. Moreover, the extension to the more realistic case Nf = 2 + 1 + 1 seems to 
be straightforward (see For studying the thermodynamic state equation, i.e. for 

computing the pressure or the energy density, as well as the mesonic screening masses 
as a function of the temperature we do not expect serious complications related to the 
presence of the twisted mass term. However, in as far the explicit flavor symmetry 
breaking will affect the thermodynamic behavior remains to be seen. 



Acknowledgements: We thank Mike Creutz, Steve Sharpe, Andre Sternbeck and 
Carsten Urbach for many useful discussions and again Carsten Urbach for his supply, 
maintenance and support of the HMC code. We also thank David Schulze for conducting 
the finite size check on L = 24. Most simulations have been done on the apeNEXT 
computers in Rome and DESY-Zeuthen. This work has been supported in part by 
the DFG Sonderforschungsbereich/Transregio SFB/TR9. O. P. and L. Z. acknowledge 
support by the DFG project PH 158/3-1. E.-M. I. was supported by DFG under contract 
FOR 465 / Mu932/2 (Forschergruppe Gitter-Hadronen-Phanomenologie). 



17 



Appendix: Simulation parameters 



In Table [T] we give a survey of all parameter sets for the full scan of the phase diagram, 
ordered according to (3. The explored k range is given as well as the values of the twisted 
mass parameter /Jo- The last four entries show the (3 range for scans performed at fixed 

(K,/Jo)- 

In Table [2] we list the locations of the identified thermal transitions K t and the 
observables, from which they were extracted. For the large (3 values the estimate is 
based on the distributions of real and imaginary parts of the Polyakov loop. If two 
values of n t are given for the same (3 and fio, this indicates that a sequence of two 
transitions (interpreted as entrance into and exit out of the cone) could be identified. 
For (3 = 3.75, 3.8 and 3.9 the critical hopping parameter for T = (found by the ETM 
Collaboration [24] at /xq = 0) is given for comparison. It should be between the lower 
and upper K t for /i ^ 0. 
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p 


K 




1.8 


0.18 < k < 0.25 


0.010/(2k) 


1.8 


0.18 < k < 0.25 


0.005/(2k) 


1.8 


0.18 < k < 0.25 


0.003/(2k) 


3.0 


0.18 < k < 0.22 


0.010/(2k) 


3.0 


0.203 < k < 0.204 


0.008/(2k) 


3.0 


0.203 < k < 0.204 


0.005/(2k) 


3.0 


0.203 < k < 0.204 


0.003/(2k) 


3.0 


0.2035 


0.0019/(2k) 


3.4 


0.16 < k < 0.24 


0.010/(2k) 


3.4 


0.16 < k < 0.20 


0.005/(2k) 


3.4 


0.16 < k < 0.1875 


0.0025/(2k) 


3.45 


0.179 < k < 0.1835 


0.0025/(2k) 


3.65 


0.169 < k < 0.176 


0.0068 


3.75 


0.1640 < k < 0.1689 


0.005 


3.775 


0.1640 < k < 0.1660 


0.005 


3.8 


0.1627 < k < 0.1675 


0.005 


3.9 


0.1550 < k < 0.1625 


0.005 


3.75 


0.1645 < k < 0.1685 


0.007 


3.75 


0.1649 < k < 0.1677 


0.008 


3.75 


0.1655 < k < 0.1678 


0.010 


3.75 


0.1655 < k < 0.1680 


0.012 


3.75 


0.1655 < k < 0.1670 


0.014 


3.75 


0.1650 < k < 0.1680 


0.025 


3.75 


0.1655 < k < 0.1675 


0.035 


3.4 


0.20 < k < 0.26 





3.5 


0.20 < k < 0.26 







U.zl < K < U.ZO 


U.UUbo 


3.45 


0.21 < k < 0.25 


0.0068 


3.65 


0.18 < k < 0.24 


0.0068 


3.75 


0.178 < k < 0.22 


0.005 


3.75 < (3 < 4.2 


0.17 


0.005 


3.7 < 13 < 5.0 


0.0001 


0.005 


3.9 < p < 4.6 


0.05 


0.005 


4.0 < (3 < 4.5 


0.1 


0.005 



Table 1: Parameter sets for the simulations. 
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p 




th 


k c (T = 0) 


observable 


3.75 


0.1656(8) 


0.005 


0.1660 


x(P) 
kl 2 
x(P) 

x(Re(L)) 
|vr| 2 

x(P) 


3.75 
3.75 
3.75 
3.75 
3.75 


0.1677(24) 
0.1657(5) 
0.1675(4) 
0.16574(24) 
0.1658(4) 


0.005 
0.007 
0.007 
0.008 
0.010 




3.775 


0.1645(1) 


0.005 




x{P) 


3.8 


0.16361(4) 


0.005 


0.164111 


|vr| 2 


3.8 


0.16621(5) 


0.005 




x(Re(L)) 


3.9 


0.1597(3) 


0.005 


0.160856 


x(Re(L)) 


3.95(10) 


0.17 


0.005 




x(Re(L)) 


4.50(5) 


0.0001 









4.45(5) 


0.05 









4.43(5) 


0.1 










Table 2: Parameter values for identified transition points. 
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